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We have developed a theoretical model and a numerical code for stationary rotating superfluid 
neutron stars in full general relativity. The underlying two-fluid model is based on Carter's covariant 
multi-fluid hydrodynamic formalism. The two fluids, representing the superfluid neutrons on one 
hand, and the protons and electrons on the other, are restricted to uniform rotation around a 
common axis, but are allowed to have different rotation rates. We have performed extensive tests of 
the numerical code, including quantitative comparisons to previous approximative results for these 
models. The results presented here are the first "exact" calculations of such models in the sense 
that no approximations (other than that inherent in a discretized numerical treatment) are used. 
Using this code we reconfirm the existence of prolate-oblate shaped configurations. We studied the 
dependency of the Kepler rotation limit and of the mass-density relation on the relative rotation 
rate. We further demonstrate how one can simulate a (albeit fluid) neutron-star "crust" by letting 
one fluid extend further outwards than the other, which results in interesting cases where the Kepler 
limit is actually determined by the outermost but slower fluid. 



INTRODUCTION 



The aim of this work is to calculate fully relativistic 
stationary models of superfluid neutron stars including 
all non-dissipative couplings between the two fluids in- 
duced by the equation of state (EOS), in particular the 
entrainment effect. In addition to studying the stationary 
properties of relativistic superfluid neutron stars, these 
models can serve as the unperturbed initial state in a dy- 
namical study of neutron star oscillations, neutron star 
collapse to a black hole, or as a starting point in studying 
pulsar glitch-models. 

Neutron stars are fascinating astrophysical objects: on 
one hand they represent a formidable "laboratory" of fun- 
damental physics, as the composition and equation of 
state of their inner core still lies beyond the reach of ex- 
perimental and theoretical physics. On the other hand, 
the advent of increasingly sensitive gravitational wave 
detectors promises to open a new observational window 
on neutron stars, which will allow us to gain new in- 
sights into these still rather poorly understood objects. 
Gravitational wave astronomy could represent the first 
opportunity to observe neutron star oscillations, provid- 
ing a new view on their inner dynamics. Considering 
the success of classical terrestrial seismology and astero- 
seismology of the sun and of main-sequence stars, one 
could expect this to result in substantial progress in our 
understanding of the dynamics and composition of neu- 
tron stars. 



tional waves 1 will give valuable complementary informa- 
tion about their rotational behavior, which is currently 
only observable via their electromagnetic pulses. 

Most theoretical studies of neutron star dynamics have 
relied on rather simplistic single-fluid models. In this 
work we attempt a more realistic description of neutron 
stars by taking their superfluidity into account via the 
use of a two-fluid model. Neutrons and protons in neu- 
tron stars are predicted to be superfluid (e.g. see [2, 3]), 
and this feature forms a fundamental ingredient in the 
current (albeit rudimentary) understanding of the glitch 
phenomenon observed in pulsars (e.g. see [4-6]). Due to 
the superfluidity and therefore lack of viscosity of the 
neutrons in the crust and in the outer core, they can flow 
freely through the other components. The remaining con- 
stituents (i.e. crust-nuclei, electrons, muons and protons) 
are assumed to be "locked" together on short timescales 
by viscosity and the magnetic field. Thereby they form 
another fluid, which in the following will be referred to as 
"protons" for simplicity. These assumptions characterize 
the so-called two-fluid model of neutron stars. These two 
fluids are strongly coupled by the strong nuclear force 
acting between protons and neutrons, and therefore a 
hydrodynamic two-fluid framework incorporating these 
couplings is required for their description. This frame- 
work will be presented in the next section. Recently it 
was pointed out that such a two-fluid system can be sub- 
ject to a two-stream instability if the relative velocity 



Additionally, observing quasi-permanent quadrupolar 
deformations ("mountains") on neutron stars via gravita- 



1 Note that this search has already begun, see [1] for a discussion 
and first results 
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of the two fluids exceeds a critical velocity [7]. This 
could therefore be relevant in neutron stars and might 
be related to the glitch phenomenon [6], which provides 
another motivation for studying the properties of such 
two-fluid systems. 

In this paper we study the stationary structure of such 
two-fluid models, in which the two fluids are restricted 
to uniform rotation around a common axis, but allow- 
ing for two different rotation rates. This neutron star 
model was first studied quantitatively by Prix [8] in the 
Newtonian context using a generalized Chandrasekhar- 
Milne slow-rotation approximation, and neglecting the 
direct interactions between the two fluids. Andersson 
and Comer [9] used Hartle's variant of the slow-rotation 
approximation to study this model in general relativity. 
Prix et al. [10] further extended the Newtonian study to 
fully include all (non-dissipative) couplings via entrain- 
ment and the nuclear "symmetry-energy", and they found 
an analytic solution for a subclass of two-fluid equations 
of state (which generalizes the P oc p 2 -type polytropes). 
More recently, Yoshida and Eriguchi [11] have devised 
an alternative approach in the Newtonian case, by treat- 
ing only the relative rotation between the two fluids as 
small, while allowing for fast rotation of the neutron star 
as a whole. Furthermore, Comer [12] has recently used 
the relativistic slow-rotation approximation to study the 
properties of the first available fully relativistic two-fluid 
EOS incorporating entrainment, which was derived by 
Comer and Joynt [13]. 

Here we present a generally relativistic numerical code 
for solving the full two-fluid model without approxima- 
tions. A preliminary progress- report on the development 
of this code, and some early results were presented in 
[14]- 

While our model and code allow in principle for any 
given two-fluid equation of state (EOS), for the sake of 
simplicity and a better numerical convergence we restrict 
ourselves in this paper to the use of a (rather general) 
class of two-fluid "polytropes". This choice is also moti- 
vated by the lack of a useful two-fluid neutron star equa- 
tion of state in the literature, especially concerning the 
aspect of entrainment. Even though Comer and Joynt 
[13] have a fully relativistic model that includes entrain- 
ment, it has not yet been developed to the point that 
it will produce a tabular equation of state that could be 
used in our code. We expect the qualitative features of 
our model to be well represented by the analytic EOS 
used in this work. 

The plan of this paper is as follows: In section II we in- 
troduce the formalism and notation of covariant two-fluid 
hydrodynamics. In section III we discuss the specializa- 
tion to an axisymmetric and stationary system, and we 
introduce the 3 + 1 framework for Einstein's equations. 
In section IV we describe the numerical procedure for 
solving the resulting elliptical system of equations. The 
tests performed on the numerical code are discussed in 
section V, and our numerical results are presented in sec- 
tion VI. A discussion of this work is given in section VII. 



In appendix A we derive a new analytic Newtonian slow- 
rotation solution, which is used for comparison to our 
numerical results. 

II. CANONICAL TWO-FLUID 
HYDRODYNAMICS 

The general relativistic framework for describing a cou- 
pled two-fluid system has been developed by Carter, Lan- 
glois and coworkers [15-18], based on an elegant varia- 
tional principle. The same relativistic two-fluid model 
was used by Andersson and Comer [9] in their slow- 
rotation description of superfluid neutron stars. 

We consider a system consisting of two fluids, namely 
neutrons and "protons", which we label by n and p re- 
spectively. The kinematics of the two fluids is described 
by the two conserved particle 4-currents n" and up", i.e. 

V a < = , and V a n* = . (1) 

The dynamics of the system is governed by a Lagrangian 
density of the form A(n",rip). Due to the requirement 
of covariance, the scalar density A can only depend on 
scalars, and we can form exactly three independent scalar 
combinations out of n" and n p , for example 

n l = --^dapn^nP, (2) 

x 2 = --^g a/3 n2nP, 

where g a p is the spacetime metric, so the Lagrangian 
density can be written as 

A«,n«) = -£(n» 2 ), (3) 

where £ is a thermodynamic potential representing the 
total energy density of the two-fluid system, or "equation 
of state". Introducing the 4-velocities it", u p of the two 
fluids, which satisfy the normalization conditions 

9af3 K u n = . and 9a/3 Up = ~C 2 , (4) 

the particle 4-currents can be written as 

n™ = n n u" , and n p = n p up" , (5) 

in terms of the neutron- and proton densities n n and 
n p respectively. Variation of the Lagrangian density (3) 
with respect to the particle currents n™ and rip 1 defines 
the conjugate momenta and namely 

dA=pldnZ+ P Zdn p '. (6) 

Due to the covariance constraint (3) we can further ex- 
press the conjugate momenta in terms of the currents 

as 

p™ = /C nn n na + /C np n pa , 

Pi = lC^n na +JC™n pa , (7) 
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where the symmetric "entrainment matrix" JC XY is given 
by the partial derivatives of £(n 2 ,n 2 , x 2 ), namely 2 



2 d£ 

c 2 dn 2 



K pp 



2 <9£ 
c 2 <9n 2 ' 



/C np 



1 d£ 

c 2 dx 2 



(8) 



The equations of motion for the two fluids can be ob- 
tained from the variational principle developed by Carter 
[19]. In the absence of direct dissipative forces acting be- 
tween the two fluids (e.g. see Langlois et al. [18]), the 
equations of motion will then be found as 3 



<V [Q ^=0, and n«V [Q p p =0 



(9) 



The energy-momentum tensor T al3 , which is derived 
from the variational principle too, has the form 



(10) 



If the equations of motion (1) and (9) are satisfied, the 
stress-energy tensor automatically satisfies V a T al3 = 0, 
which is a Noether-type identity of the variational prin- 
ciple. The generalised pressure '5 of the two-fluid system 
is defined by the thermodynamical identity 



We can now equivalently consider the equation of state 
£ as a function of the form £(n n , n p , A 2 ), for which the 
first law of thermodynamics reads as 

d£ = n a dn n + fjPdn p + a dA 2 , (15) 

closely analogous to the Newtonian formulation [10]. The 
conjugate quantities defined in this equation are the en- 
trainment a and the neutron- and proton chemical po- 
tentials ^t n and /i p (sometimes also referred to as specific 
enthalpies, which is equivalent in the zero-temperature 
case). It is often useful to characterize the entrainment 
by the dimensionless entrainment numbers ex , which we 
define as 



2a 



m x nx 



(16) 



where m x is the particle rest-mass of the respective fluid, 
and the fluid- index is X = n, p (no summation over X) . 

The conjugate variables of (15) can be expressed in 
terms of the kinematic scalars and the entrainment ma- 
trix JC XY as 



(/C nn n 2 +/C np .T 2 ) =-<pS, 



£ + * 



-<p n a 



<P p a , 



(11) 



so '5 can be considered the Legendre-transform of the 
Lagrangian density A. Using the entrainment relation 
(7), we can rewrite this as 



£ + # 



= /C nn n 2 + 2/C np x 2 + JC PP n 2 . 



(12) 



Instead of x 2 defined in (2) , we will use a physically more 
intuitive quantity as the third independent scalar, namely 
the "relative speed" A. We define the relative speed A 
as the norm of the neutron velocity as seen in the 
frame of the protons u p , or vice versa. The corresponding 
relative Lorentz factor Ta is therefore given by 



1 x 2 ( 

p n n n p V 



-1/2 



(13) 



tr r n n n p \ c 

and the relative speed A is expressible in terms of x as 

x" 



_ ( n n n p y 
\ x 2 J 



(14) 



In the case of co-moving fluids (i.e. — u p ), we see 
from (2) that x 2 = n D n p , and so A = as expected. 



M p = — (!C pp n 2 +/C np a 



a = 



c 

i/C np n n n p ri . 



u pl J a 



(17) 



Using (14), the inverse relations can be obtained as 



/C nn = 
/C pp = 
/C np = 



n n <r 



„2p2 

2a 

„2p2 
"p 1 A 



(18) 



2a 



r 3 ' 

1 A 



which reduces exactly to the corresponding relations in 
the Newtonian limit [10], where Ta 1 and ^ — ► 
m x c 2 +p x . In terms of these quantities, the generalised 
pressure ^ (11) can also be written as 



-£ + n n fi n + n p /j, p , 



(19) 



which is the generalization of the thermodynamic Gibbs- 
Duhem relation to a two-fluid system. 



III. STATIONARY AXISYMMETRIC 
CONFIGURATIONS 



2 The corresponding notation in Andersson and Comer [9] is ra™ — » 
n a , n£ -» p a , < -> u a , u« -> v a , K n P -> A, K nn B, and 
Cpp -» C. 

3 the square brackets denote averaged index anti-symmetrization, 
i-e- 2i) |tli i,] = v ab - v ba . 



A. The metric 

Here and in the following we choose units such that 
G = c = 1 for simplicity. We consider spacetimes that are 
stationary, axisymmetric, and asymptotically flat. The 
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symmetries of stationarity and axisymmetry are associ- 
ated with the existence of two Killing vector fields, one 
timelike at spatial infinity, t a , and one spacelike every- 
where and with closed orbits, (p a . 

It was shown by Carter [20] that under these assump- 
tions the Killing vectors commute, and one can choose 
an adapted coordinate system (t, x , x 2 , ip), such that 
t a d a = d/dt and p a d a — d/dp, i.e. 

t a = (1,0,0,0), and p a = (0,0,0,1). (20) 

We choose the remaining two coordinates to be of spher- 
ical type, i.e. x 1 = r, x 2 = 9, and following Gourgoulhon 
et al. [21], we fix the gauge to be of maximal-slicing quasi- 
isotropic type (MSQI), for which the line element reads 
as 

ds 2 = g al3 dx a dx fi = —(N 2 - N^N^) dt 2 - 2N V dip dt 
+A 2 (dr 2 + r 2 d9 2 ) + B 2 r 2 sin 2 <9 dp 2 , (21) 

where the functions N, N v , A and B depend on r and 9 
only, and N v = g vv N v . 

B. Fluid dynamics 

We assume the flow of the two fluids to be purely axial 
(i.e. no convective meridional currents), so we can write 
the unit 4-velocities of the two fluids as 

< = and = (22) 

where the helical vectors and £p are expressible in 
terms of the Killing vectors t a and p a as 

C=t a + n n p a , and Cp=t a + %p a , (23) 

where the two rotation rates O n and f2 p are scalar func- 
tions, which can only depend on r and 9. 

Using Cartan's formula for the Lie derivative of a 1- 
form p/3 with respect to a vector-field £ a , namely 

Ct Pa = 2^Vpp a] + V Q (f P0 ) , (24) 
we can rewrite the equations of motion (9) as 

A*F£-Va(Cj?p£) =0. (25) 

Linearity of the Lie derivative together with (23) and (24) 
allows us to rewrite this as 

CtpZ + nxCvP^ + pPpXpVaflx-Va^pXp) =0. (26) 

Stationarity and axisymmetry imply that the first two 
terms vanish, and so the equations of motion for neutrons 
and protons are reduced to 

p£V a fi„ = V Q {C P n p) , <V Q fi p = V Q (C^) .(27) 



In the general case of differential rotation, the integrabil- 
ity condition of these equations are therefore 

p£=j$(fi„), and pP=pP(tt p ), (28) 

and the first integrals of motion are obtained as 

p%(Q.')dQ,' = const 11 , (29) 

V\ + n pP l - pP (n')dfi' - const? . (30) 

In the special case of uniform rotation, i.e. V a flx = 0, 
these first integrals reduce to 

p n t + ft n p£ = const 11 , and p\ + Q p pP = const p , (31) 

which are equivalent to the expressions obtained by An- 
dersson and Comer [9]. We can further express these 
first integrals in terms of the chemical potentials [i n , p p 
of (17), namely 

Pt + ^P% = CpI = — yA 1 " = const 11 , (32) 

u n 

p? + n v pl = C P >S = = constP . (33) 

C. The 3 + 1 decomposition 

We introduce the vector n" as the unit normal to the 
spacelike hypersurfaces S t defined by t = const., namely 

n a = -NV a t, (34) 

which defines the so-called Eulerian observers Oo follow- 
ing Smarr and York [22]. The induced metric h a p on the 
spacelike hypersurfaces S t is given by the projection 

h a /3 = 9afi + n a np . (35) 

The corresponding 3 + 1 decomposition of the energy- 
momentum tensor T a/3 reads as 4 

T a(3 = S a(3 + 2n {a J r3) + En a n p , (36) 

where 

E = n^Tapxi , J a = -h^Tj^n^ , S a p = h^T^^hp , 

(37) 

which can be interpreted as the energy, momentum and 
stress tensor as measured by the Eulerian observers. In 
the MSQI gauge (21), we can explicitly express these 
quantities as 

E = N 2 T U , Ji = NT* , S) = T) - N l T] , (38) 



4 Round brackets denote averaged symmetrization, i.e. 2u( a j,) = 

Vab + v ba- 
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where Latin indices i, j denote the space-components 
1,2,3. The Einstein equations in this formulation result 
in a set of four elliptic equations for the metric potentials 
(see [23] and [21] for details), namely 

A 3 v = 4irA 2 (E + SI) + A 2 K l0 K^ - dvd{v + (3), (39) 
A 3 N V = -WirNA 2 ^ - rsm6dN v d(3(3 - v) , (40) 
A 2 [(NB - l)rsin0] = 8nNA 2 Brsm6 (S r r + S e g ) , (41) 

A 2 {v + a) = 8ttA 2 S% + Z -A 2 K l3 K^ - {dvf , (42) 

where we defined N v = rsinON^ and J v = rsin6> J v . 
A3 and A2 are the flat three- and two-dimensional 
Laplace operators, whereas A3 = A3 — (r 2 sin 2 O)^ 1 . We 
further used the notation 



+ -dad/3 
8 r (a - (3/2) + 









dV + 











1 



r tan# 



d e (a - 0/2) 



. (49) 



Both of these virial theorems (48) and (49) can be written 
as the sum of an integral over a "material" term / ma t 

(the 

first term in (48) and (49) respectively), and an integral 
over pure field-quantities /fields (the remaining terms). 
Therefore it will be convenient to consider the following 
dimensionless quantity to numerically characterize the 
respective virial violations: 



GRV 



elds 



(50) 



v = \nN . 



a 



In A, (3 = \nB, 



(43) 



and we define dad/3 as the flat-space scalar product of 
two gradients, i.e. 



dad/3 = d r ad r /3 + \deadef3 . 



(44) 



The only non-zero components of the extrinsic curvature 
Kij in our spherical coordinate basis are given by 

K rv = -g^ d r N* , and K dip = -f^ d e N* . (45) 

We note that the gravitational mass Ai , which is defined 
as the (negative) coefficient of the term 1 jr in an asymp- 
totic expansion of the "gravitational potential" log TV, can 
be expressed explicitly (see [23]) as 



M 



I 



A 2 B N(E + S t i ) + 2B 2 N ip J v r 2 sin 9 dr d9 dip . 

(46) 

Here and in the following we will use M. to denote the 
gravitational mass, while M will stand for the baryon 
mass. The total angular momentum J is given by 



J 



-I 



A 2 B 3 r sin9 J v 



r 2 sin 9 drdO dip. (47) 



The 2D- and 3D virial identities, which have been derived 
by Bonazzola and Gourgoulhon [24, 25], can serve as a 
useful check of consistency and precision of the numerical 
results. The 2D virial identity (referred to as GRV2), 
which derives from the Poisson-equation (42), has the 
form 



/ 



8ttA 2 S% + -A 2 KijK lj - (dvf 



rdrd6 = 0, (48) 



while the 3D virial identity (GRV3), which reduces to 
the usual virial theorem in the Newtonian limit, can be 
written as 



/ 



AwA 2 BS l i dV + / B 



4 3 



{dvf 



D. The matter sources 

Let us write T n and T p for the two Lorentz factors 
linking the Eulerian observers Oq to the co-moving fluid 
observers O n (defined by w") and O p (defined by up 1 ), 
namely 

T n = — n a u" = Nu^ , and T p = — n a Up = Nu P . 

(51) 

The "physical" fluid velocities U n and U p of the two flu- 
ids 5 in the ip direction, as measured by Co, are given 
by 

U n = YTVaK . and U P = TTVaU p , (52) 

where cp a is the spatial unit vector in the ip direction, i.e. 
1 



p a , such that h af3 <p> a p p = 1 . (53) 



Using (22) and (51), we obtain 



U n 



(O n - iV") , U p 



(O p - N") , (54) 



N y " ' ' " N 
and the Lorentz factors can be expressed equivalently as 

-1/2 



Tn = (i-u 2 n y 



and r r 



1-UlY 112 . (55) 



The "crossed" scalar x 2 , defined in (2), can be expressed 
in terms of the respective scalar particle densities n n , n p 
and the 3-velocities U n and U p , as 



U n U p 



^(1-U 2 )(1^U 2 )' 



(56) 



6 In Andersson and Comer [9] these were denoted — u)n and — lu p 
respectively. 
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and using (14), we can write the relative velocity A as 



(Un ~ U P ) 2 

(1 - C/ n C/ p ) 2 



(57) 



Using expressions (32), (33) and (51), the first integrals 
can be cast into the form 



N n n , N n 

— n = const , and 

1 n 1 p 



const p . 



(58) 



In closer analogy with [21, 23], we can alternatively write 
these first integrals as 



H n + v-lnT n = C n . 
H p + u-]nT v 



■'p ' 



where we introduced the abbreviations 



H n = ln 



(—) 

\m n J 



and H n 



ln(^ 

\m p 



(59) 
(60) 



(61) 




The components of the 3 + 1 decomposition (36) of the 
energy-momentum tensor (10) are explicitly found as 

E = -* + {T 2 n }C nn nl + T 2 p K, pp nl 

+2r n r p /C np n n n p ) , (62) 

= rl>c nn nlu n + rl>c pp nlu p 

+T n T p IC np n n n p (U n + U p ), (63) 
- $ = (64) 
= * + (T 2 n lC™n 2 n U* + T 2 p IC pp nlU 2 

+2T n T p IC np n D n p U n U p ) , (65) 

One can check the consistency of this result with the 
single fluid case of [23] , by considering the special case of 
both fluids moving together. 



IV. NUMERICAL PROCEDURE 
A. Iteration scheme 

The numerical solution of the stationary axisymmetric 
configurations described in the previous sections pro- 
ceeds in a very similar manner to the single-fluid case, 
which is described in more detail in [21, 23]. The central 
iteration scheme is nearly identical: 

Initialization: Start from a simple "guess" for a spher- 
ically symmetric matter distribution and n p °' 
of the two fluids, and use a flat metric. 

Step 1: Calculate the matter source-terms E, J v and 5j 
from (62)- (65). 

Step 2: Solve the equations (39)-(42) for the corre- 
sponding metric using the pseudo-spectral elliptic 
solver in LORENE, the numerical relativity pack- 
age used here[26]. 



Step 3: Use the first integrals (58) to obtain the chemi- 
cal potentials /j, n and n p . 

Step 4: Calculate the new density fields n n and n p by 
inverting the relations (17) for the given equation 
of state. 

Step 5: Continue at Step 1 until the desired convergence 
is achieved. 

In general we are using three spherical-type numerical 
domains to cover the hypersurface : the innermost do- 
main covers the whole star. An intermediate domain is 
used for the vicinity of the star to about twice the stellar 
radius, and an outer domain covers the space out to infin- 
ity, using a compactification of the type u = 1/r (see [23] 
for details). For the inner domain we have the choice of 
either using a simple spherical grid containing the whole 
star, or we can use an adaptive-grid algorithm, in or- 
der to adapt the domain-boundary to the stellar surface. 
Contrary to the single-fluid case (e.g. see [21]), however, 
the adaptive-grid approach is much less effective in in- 
creasing precision and convergence. The reason for this 
is simply that only one of the two fluid-surfaces can be 
matched up with a domain boundary, and therefore the 
(weak) Gibbs phenomenon due to the inner fluid surface 
(representing at least a discontinuity in the derivative) is 
not completely avoidable. 

Another important difference in the case of two-fluids 
compared to single-fluid stars is the way we determine 
the location of the fluid surfaces. In a single fluid star, 
the surface can always be found by the vanishing of the 
pressure, which usually translates into a simple condition 
in terms of the vanishing of the chemical potential /i. In 
the two-fluid case, however, this is not generally possible 
(especially for the "inner" fluid), due to the coupling of 
the fluids. Therefore we need to define the fluid surfaces 
directly in terms of the vanishing of the respective density 
fields. Contrary to the chemical potential, the density 
can have a vanishing or diverging gradient at the surface, 
and a precise determination of the surface can therefore 
be numerically difficult. 6 

A related numerical problem specific to two-fluid con- 
figurations appears when the surfaces of the two fluids are 
very close to each other. In this case, the 1-fluid region in 
between the two surfaces will be poorly resolved by the 
grid covering the star, and therefore the determination 
of the outer surface will have a poor numerical precision. 
As will be seen later, this problem can be cured to some 
extent by adding another domain, which covers just a 
thin shell below and up to the outer fluid surface. In this 
case one observes a drastic improvement in the precision 



We note that Yoshida and Eriguchi [11] chose to avoid this dif- 
ficulty by defining the "fluid surfaces" by the vanishing of the 
respective chemical potentials. These "surfaces", however, do 
generally not coincide with the surfaces of vanishing density (con- 
trary to the single-fluid case), as can be seen from (17). 
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of finding the outer surface, which can be quantified by 
comparison with the analytic slow-rotation solution. 

We note that this numerical code can be used equally 
well for Newtonian configurations, simply by replacing 
the matter-sources by their Newtonian limits, and forc- 
ing the spatial metric to be flat. The central iteration 
scheme remains unchanged, and we can relate the lapse 
N to the Newtonian gravitational potential, namely by 
the relation v — In AT = Q/c 2 , where $ is the Newto- 
nian gravitational potential. The Newtonian limit of the 
matter source-term in Eq. (39) is 

^l=p + 0(c- 2 ), (66) 

where p is the total (rest-)mass-density, so that this com- 
ponent of the Einstein equations reduces to the New- 
tonian Poisson equation, while the remaining Einstein 
equations (40)-(42) become trivial in this limit. In a sim- 
ilar manner, the first integrals are seen to reduce exactly 
to their Newtonian counterparts [10]. 

The parameters of the numerical scheme that will be 
used for the rest of the paper are the following: the re- 
quired convergence of the iteration scheme is 10 -10 , and 
we use 17 points in the 9 direction, and 33 grid-points 
in the radial direction in the innermost domain (contain- 
ing the star), 33 radial points in the intermediate domain 
and 17 radial points in the compactified outer domain. 

B. The polytropic 2-fluid equation of state 

The numerical scheme described in the previous section 
can be used for any invertible 2-fluid equation of state 
(EOS). The current implementation of our code, however, 
only covers a "polytropic" subclass of 2-fluid EOS, which 
generalizes the types of EOS used in previous studies, 
e.g. [9-11, 27], and which has the general form 

£ = pc 2 + ^ n n 71 + ^k p ^ 2 + K np < 3 np 4 

+K A n?> p ' 6 A 2 , (67) 

where p = m n n n + m p n p . As discussed in the intro- 
duction, we expect this polytropic EOS-class to be quite 
general, and to allow one to study the qualitative fea- 
tures of a broad range of different superfluid neutron star 
models. For example, general features of the Kepler limit 
(cf. Fig. 5) are seen to be in qualitative agreement with 
the mean field results of Comer [12]. 

The two fluids in (67) are coupled via a "symmetry 
energy"-type term proportional to K np and an entrain- 
ment term proportional to k a A 2 . The resulting expres- 
sions for the chemical potentials and the entrainment a 
are directly obtainable from (15). 

In general this class of 2-fluid EOS requires a numer- 
ical inversion in the iteration scheme described in sec- 
tion IV A, in order to obtain the densities n n ,n p from 
the chemical potentials p n , pP at a given relative speed 



A. For testing purposes and for comparison to the New- 
tonian and relativistic slow-rotation results, we will in 
the following be mostly interested in a further subclass 
of the above EOS, namely the special 2-fluid polytropes 
described by 

£ = P c 2 + ^K n nl + ^K p nl + K np n n n p + K & n n n p A 2 , (68) 

which are a 2-fluid generalization of the 1-fluid polytrope 
P oc n 2 . This special EOS class still exhibits all the 
coupling-types (entrainment + symmetry energy) of the 
general EOS, but allows an analytic inversion, namely 

p n -m n c 2 = K n n n + (k„ p + k a A 2 ) n p , (69) 
p p -m p c 2 = K p n p + (n np + k a A 2 ) n n , (70) 

and the entrainment is found as 

a = K A n n n p . (71) 

The generalized pressure in (19) is expressible as 

* = 2 KnB - 2 + 2 KpTl P + Kn P nnTl P + K A n n«- P A 2 . (72) 

Contrary to the two-fluid EOS used in the Newtonian 
slow-rotation study [10], which exhibits the somewhat 
unphysical feature of constant entrainment numbers, as 
discussed in appendix A, this EOS results in a much 
more physical behavior of the entrainment. Namely, us- 
ing (16), we find 

£ n = n p , and e p = — ^ n n , (73) 
in m p 

which ensures that the entrainment effect automatically 
vanishes when one of the two fluid-densities goes to zero. 
Such a linear behavior of entrainment also happens to be 
in quite good qualitative agreement with the theoretical 
predictions of nuclear physics, e.g. see 7 [2, 3, 28] 

Using the method developed in [10] for the EOS (68), 
we can find an analytic solution in the Newtonian slow- 
rotation approach, which is presented in appendix A. 
This allows us to run extensive tests by comparing our 
numerical code to the analytic solution in the Newtonian 
case. The results of this comparison are presented in 
section V C. 



V. TESTS OF THE NUMERICAL CODE 

A. Comparison to 1-fluid results 

As a first consistency check we use the two-fluid code 
for strictly co-rotating configurations with a common 



7 These references give the neutron and proton effective masses 
m"*, which are related to the entrainment via ex = (rre* — 
ir^'l/m 1 , see [10] for details. 
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outer surface, and compare the results to those of the 
well-tested single-fluid code [21, 23]. For this purpose we 
study a stellar sequence of fixed central density and vary 
the rotation rate. We define the "natural scale" of the 
rotation-rate as 



O = v/47rGp(0), 



(74) 



where p(0) is the central rest-mass density, i.e. p(0) = 
m n n n (0) + m p n p (0). The Kepler rotation rate Ok is typ- 
ically found at about Ok ~ 0.1 Oo for the configurations 
considered here. The results of the comparison with the 
single-fluid case are shown in Fig. 1. Here we plot the 
relative differences, defined as 



Q 11 



(75) 



of a global quantity Q in the two-fluid case (Q ) and in 
the single-fluid case (Q lf ). The first column, figures 1 (a) 
and (c), shows the comparison of 1-fluid and 2-fluid re- 
sults using a fixed spherical grid for the inner domain 
in the two-fluid case. The single-fluid code on the other 
hand always uses an adaptive grid for the stellar surface. 
We notice that towards higher rotation rates the relative 
errors increase. These errors can be entirely ascribed to 
the lack of grid adaption in the two-fluid case: by using 
an adaptive grid for the stellar domain also in the two- 
fluid case, we find a consistent agreement of better than 
1CP 9 , as can be seen in the second column in figure 1 (b) 
and (d). We note, however, that this improvement of 
using an adaptive grid is restricted to cases where the 
two-fluids share a common outer surface, while it is of 
much less use in the general two-fluid case as mentioned 
earlier. We can therefore conclude that the two-fluid code 
reproduces results consistent with the single-fluid code in 
cases where the two fluids co-rotate. 



B. Virial theorem violation 

In the next step we consider the more general case 
where the two fluids have different rotation rates. We fix 
the relative rotation rate, defined as 



K = ~ o7~ 



(76) 



to be 1Z — 1.51 and vary O n . As mentioned before, in 
these general situations an adaptive grid does not sub- 
stantially improve the precision and has therefore not 
been used. We consider the internal consistency check 
provided by the virial identities GRV2 and GRV3 de- 
fined in (50), for which the result is shown in Fig. 2. We 
note that in the case (a) , where one inner domain is used 
to cover the star, even at low rotation rates the result 
falls somewhat short of the convergence-goal of 10~ 10 in 
the iteration scheme. This lack of precision at small ro- 
tation rates can be understood as follows: due to the 



difference in rotation rates, the two fluids do not share 
a common outer surface, and there will necessarily be a 
1-fluid region close to the outer surface. However, this 1- 
fluid region will be very thin compared to the dimensions 
of the star, and will therefore be poorly resolved in terms 
of the numerical grid. We can improve this by choos- 
ing a second domain to cover just a thin layer (of about 
1% of the radius) below the outer surface, resolved by 
another 33 radial grid-points. The effect of this "trick" 
is rather impressive and can be seen in Fig. 2, for the 
case (b). While this gain in precision is not very impor- 
tant by itself, it underlines the consistency of the results 
and shows that the source of these errors is understood. 
The decrease in precision when approaching the Kepler 
rotation can be ascribed to the appearance of cusps at 
the equator (see Fig. 6) and therefore the presence of the 
Gibbs phenomenon. Nevertheless, one should note that 
this phenomenon happens also in the one-fluid case, the 
precision of the code at the Kepler limit being of the same 
order as here [21]. 



C. Comparison to Newtonian slow-rotation results 

We can use the analytic Newtonian solution in the 
slow-rotation approximation (derived in appendix A) for 
a systematic comparison with the numerical code run in 
"Newtonian mode" as described in section IV A. 

We denote the numerical solution of a quantity as 
Ql(Ox), and the analytic slow- rotation solution as 
Q sr (Ox), and we define the relative difference as 



oQ = 



Ql — Q s 



(77) 



where corresponds to the static solution. At fixed 
relative rotation rate TZ, the slow-rotation solution can 
be written as 



0sr(O n ) = Q<°) + 02. 



(78) 



Ideally we would like to compare only up to the O 2 com- 
ponent of the numerical solution, but obviously we do not 
know its Taylor-expansion in orders of O. Nevertheless 
the numerical solution can formally be written as 



.(On) 



(79) 



so that the relative difference (77) can be expanded as 
Q<°> 



o (4) 



(80) 



If the numerical solution agreed perfectly with the ana- 
lytic solution (up to order O 2 ), the first two terms would 
be zero, and the leading order of the difference would be 
O n . In practice, however, there will be contributions on 
all orders, and we will try to quantify these respective 
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FIG. 1: Relative differences AQ in the total baryon mass M, and the equatorial and polar radii 7? eq and 7? po i. The first 
row is the Newtonian case, while the second row shows the relativistic results. In the first column, we used a fixed spherical 
inner domain, while in the second column the inner domain-grid is adapted to the stellar surface (which is the default in the 
single-fluid case). 
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FIG. 2: Virial violations GRV2 and GRV3 for (a) 1 domain 
and (b) 2 domains covering the star. Only the Newtonian 
case is shown, as the relativistic results are very similar, flu 
denotes the Kepler rotation rate. 



errors. If in some interval of fi n one of these terms dom- 
inates in the series, then a log-log plot of oQ(fl n ) in this 



region would look like 

2/ = log(afO =loga + mlogn n , (81) 

i.e. a straight line with steepness m and an offset log a. 
Conversely, if the log-log plot of oQ contains sections of 
straight lines, we can infer the leading power in Q n and 
its coefficient. 

The neutron-star model used here is characterized by 
the following choice of EOS-parameters: 

k d = 0.02 , k p = 0.12 , k„ p = 0.01 , k a = 0.02 , (82) 

and the (fixed) central chemical potentials 

M n (0) = M p (0) = 0.2 m b c 2 , (83) 

where nit, = 1.66 x 10~ 27 kg is the baryon mass and c is 
the speed of light. The resulting neutron-star model in 
the static case has a total mass of M = 1.50 M Q (where 
Mq is the solar mass), a radius of R = 11.1 km, a central 
baryon number density n(0) = 1.04 fm~ 3 , proton fraction 
x p = n p /n = 0.083 and an entrainment value of e P (0) = 
0.38. We note that in the Newtonian case there is no 
distinction between the gravitational mass M. and the 
baryon mass M. 

In Fig. 3 we show the relative differences oQ for the 
radii R n and R p at the equator. We also plotted the 
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straight lines corresponding to a pure fi 4 and fi 2 behav- 
ior, in order to facilitate the interpretation of these re- 
sults. 

In Fig. 3 (a) we see that for small rotation rates 
(Q n /Q,Q < 10 -3 ) the error in the equatorial proton radius 
i? p (eq) reaches a "plateau" at about <~ 10 -9 , which cor- 
responds to numerical errors and the finite convergence- 
condition of the iteration scheme, while for higher rota- 
tion rates, the quartic error starts to dominate. The same 
behavior is observed for other global quantities, e.g. M n 
and M p and i?(pol)), which are not included in this plot. 
However, the neutron equatorial radius i? n (eq) (which is 
the outer radius) displays a consistent quadratic error of 
order unity! The reason for this apparent discrepancy 
is rather subtle, and stems from the somewhat different 
nature of the slow-rotation approach and the fully nu- 
merical solution. In the numerical code, when one of 
the two fluid-densities vanishes, we switch from the 2- 
fluid EOS (67) to the corresponding 1-fiuid EOS before 
we do the inversion ^ — > riy in the numerical proce- 
dure (cf . section IV A) , which is the correct physical way 
to do this. In the slow-rotation approach, however, the 
rotation rates are treated as infinitesimal, and there is 
actually no finite 1-fluid region. Therefore the EOS is 
always used in the form (67), which will be seen in the 
following to account for the difference in i? n (eq). In or- 
der to test this explanation, we have also implemented a 
"slow-rotation style" EOS-inversion in the code, in which 
we do not switch to a 1-fluid EOS when one of the two 
fluids vanishes. The result of this is shown in Fig. 3 (b). 
We see that the discrepancy of the outer radius has com- 
pletely disappeared. While this serves as an interesting 
test of consistency, this rather unphysical EOS-inversion 
will obviously not be used in the following. 

D. Comparison to relativistic slow- rotation results 

Finally, we compare our results in the fully relativistic 
case to the results obtained by using a code developed 
by Andersson and Comer [9] , which is based on the rela- 
tivistic slow-rotation approximation. 

In the relativistic case the physical "radius" will gener- 
ally be different from the coordinate-radius, and can be 
defined in various non-equivalent ways (e.g. circumfer- 
ential radius, proper radius). For an unambiguous com- 
parison we define the "radius" R as the proper distance 
of the surface from the center of the star, along a line 
of constant 9 and tp (the definition of which is consistent 
with [9]), i.e. 

R= dl= A(r)dr, (84) 
Jo Jo 

where Rq is the coordinate-radius of the surface. Another 
quantity specific to the relativistic case is the shift- vector 
N 1 = (0, 0, N v ), and we will consider its 3-norm, i.e. 

\\N i \\ = ^¥m=N^^, (85) 



which is independent of the coordinate-system chosen on 
the spacelike hypersurface. 

The stellar model used in this comparison is defined 
by the EOS parameters 

n n = 0.04 , k p = 0.24 , « np = 0.02 , k a = 0.02 , (86) 

and the central chemical potentials are ^"(0) = /i p (0) = 
0.2 rribC 2 . The configurations obtained have the follow- 
ing (fixed) central values: the central baryon density is 
n(0) = 0.5776 fm~ 3 , which corresponds to 3.61 times nu- 
clear density (n nuc i = 0.16 fm -3 ). The central proton 
entrainment is £ p (0) = 0.212, and the proton fraction is 
found as x p (0) — 0.083. We fix the relative rotation to 
1Z = 0.5, i.e. the neutron superfluid is rotating 50% faster 
than the proton-electron fluid. In Table I we show the 



fi n /27T 


Hz 


100 Hz 


500 Hz 


M n [Mq] 


1.0978 (-0.02%) 


1.0998 (-0.1%) 


1.1509 (-2%) 


M p [Mq] 


0.0998 (-0.02%) 


0.0997 (-0.1%) 


0.0959 (-2%) 


M[Mq] 


1.1194 (-0.02%,) 


1.1210 (-0.04%) 


1.1644 (0.04%) 


R? [km] 


13.545 (-0.01%) 


13.570 (0.2%) 


14.260 (5%) 


RZ o1 [km] 


13.545 (-0.01%) 


13.527 (-0.1%) 


13.103 (-3%) 


R? [km] 


13.545 (-0.01%) 


13.534 (-0.1%) 


13.302 (-2%) 


R P P o1 [km] 


13.545 (-0.01%) 


13.527 (-0.1%) 


13.103 (-3%) 


N(0) 


0.700102 (le-4%) 


0.69983 (2e-3%) 


0.69267(-0.04%) 


imi(eq) 





0.00206 (0.1%) 


0.01072 (3%) 



TABLE I: Numerical results Ql and (in parentheses) relative 
differences (Ql — Qst)/Ql x 100% to the relativistic slow- 
rotation results Q sr . 

results of the comparison to the relativistic slow-rotation 
code. We observe that generally the agreement is quite 
good, and (as expected) gets worse with higher rotation 
rates. However, we note that this slow-rotation code im- 
poses an additional constraint on the radii, namely the 
two fluids are forced to share a common outer surface. 
Therefore part of the disagreement observed here does 
not actually stem from the slow-rotation approximation 
or numerical differences, but from the somewhat different 
assumptions in the model. Given these differences, the 
agreement seems very good. 

VI. NUMERICAL RESULTS 

The existence of configurations with one fluid-surface 
having a prolate shape was initially found using the New- 
tonian analytic solution[10] in the slow-rotation approx- 
imation. While this might not be very realistic astro- 
physically, it is still interesting to study this particu- 
larity of an interacting two-fluid system. We confirm 
the existence of such configurations in the fully relativis- 
tic treatment, as reported earlier by us[14]. In order 
to show this, we choose the polytropic EOS parameters 
n n = 0.016, k p = 0.16, K np = 0.008, and k a = 0.03, 
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log 10 (fi n /fio) log 10 (n n /f2o) 



FIG. 3: Relative difference oQ(Sl n ) between the numerical code in "Newtonian mode" and the slow-rotation analytic solution 
of appendix A, for the equatorial radii R n and R p . In (a) we used the normal "physical" EOS-inversion, while (b) shows the 
results when using a "slow-rotation style" EOS inversion. 



with the central chemical potentials /Lt n (0) = 0.2mf,c 2 and 
A* P (0) = 0.198mf,c 2 . This corresponds to a central proton 
fraction of x p (0) — 0.05 and a central proton-entrainment 
number of e p (0) = 0.80. In Fig. 4 we show the result- 




-5 5 

x [km] 

FIG. 4: Meridional cross-section of an oblate-prolate two-fluid 
configuration. The dotted lines represent lines of constant 
"gravitational potential" N, while the thick lines are the re- 
spective surfaces of the neutron- and proton fluids. 

ing configuration with the two fluids counter-rotating at 
n n /2w = 1000 Hz and n p /2ir = -100 Hz. We define the 



ellipticity of fluid X as 

_ Rxjcg) - Rxjpol) , . 

e x = 5-7 — \ , (87) 

Rx (eq) 

in terms of the proper radii R of (84). Using this defini- 
tion, this configuration is found to have e n = 0.137, and 
e p = —0.037, so the proton fluid has a prolate shape de- 
spite the fact that it is rotating around the z-axis. This is 
made possible by the effective interaction potential cre- 
ated by the neutron-fluid, which "squeezes" the proton- 
fluid, in this case to the point of even overcoming the 
centrifugal potential. 



EOS 










I 


0.05 


0.5 


0.025 


0.02 


II 


0.05 


0.5 


0.0 


0.0 


III 


0.05 


0.5 


-0.025 


0.02 



TABLE II: Polytropic parameters defining EOS-models I, II, 
and III 

To simplify the presentation of results, we focus in the 
following on three EOS-models, defined in table II, which 
differ only by their interaction-terms. The EOS-models I 
and III differ by the sign of the "symmetry-interaction" 
term n np , which corresponds to a value of the canon- 
ical "symmetry-energy term" (introduced in Prix et al. 
[10]) of cr = -0.5 for EOS I and a = 0.5 for EOS III. 
EOS II represents two fluids without EOS-interactions. 
If not otherwise stated, we choose the central chemical 
potentials as /j, n = /i p = 0.3 m^c 1 . In the static case we 
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FIG. 5: Kepler limit S1k as a function of relative rotation rate TZ for EOS-models I, II and III. The dashed line (NSR) represents 
the result from the analytic Newtonian slow-rotation solution (cf. Appendix A) 



obtain the results shown in Table III for these three EOS- 
models. We note that all three static configurations are 
on the stable branch of the mass-density relation, which 
can be seen in Fig. 7 for EOS-model I. We note that 





EOS I 


EOS II 


EOS III 


n c [fm 3 ] 


0.7177 


0.7697 


0.8612 


Ep(0) 


0.273 


0.0 


0.301 


i p (0) 


0.05 


0.09 


0.125 


M [M ] 


1.586 


1.532 


1.448 


M [M Q ] 


1.460 


1.409 


1.332 


R [km] 


14.37 


13.88 


13.12 



TABLE III: Results for the central baryon number density 
n c , entrainment e p , proton fraction x p , total baryon mass 
M, gravitational mass M and the proper radius R for EOS- 
models I, II and III in the static case. 

when considering rotation, the individual fluid radii will 
obviously change, but also the masses, because we only 
consider stellar sequences of fixed central density. 

Next we consider these stellar models rotating at their 
maximum rotation rate fix (called Kepler limit) for dif- 
ferent relative rotation rates TZ, as shown in Fig. 5. We 
define the Kepler-rate as the rotation rate of the 
outer fluid (which in this case also happens to be the 
faster rotating one), i.e. the protons for TZ < and the 
neutrons for TZ > 0. The rotation rate of the (slower) 
inner fluid is trivially determined by f^K and TZ. The 
dashed line shows the result from the Newtonian slow- 
rotation solution (A42). This is seen to overestimates 
the Kepler rate typically by about 15 %, except for the 
case of EOS I, where it can even underestimate the Ke- 
pler limit for TZ < 0. We see that in the fixed central- 
density sequences considered here, the local maximum 
of the Kepler rate is always attained for the co-rotating 



configuration (i.e. TZ — 0), which contrasts with the case 
of fixed-mass sequences as considered in the Newtonian 
study [10]. A similar feature of the Kepler-rate decreasing 
as TZ decreases through zero can be seen in the mean-field 
results of Comer [12]. 

In the Fig. 6, we show the fluid surfaces of the two 
fluids rotating at the Kepler-rate for two different rela- 
tive rotation rates, TZ = 0.1 and TZ = 0.01 respectively. 
We see the characteristic "cusp" appearing at the equa- 
tor of the outer fluid, which indicates the onset of mass- 
shedding if the rotation rate were to be increased any 
further. Because we fixed the central densities of these 
configurations to those of the static case, it can be seen 
from Fig. 7 that both of these configurations belong to 
the so-called "supramassive" class, i.e. they do not have a 
corresponding stable non-rotating configuration of equal 
baryon-mass. 

Fig. 7 shows the mass-density diagram for the static 
configuration of EOS I and for three Kepler configura- 
tions with different relative rotation-rates. The configu- 
rations to the right of the maximum are on the so-called 
"unstable branch", because they will be subject to un- 
stable modes under small perturbations. The configura- 
tions above the dotted line correspond to stars on the 
unstable-branch of the static curve. They have no sta- 
ble non-rotating counter-part, even if they are on the 
stable branch of the mass-curve of the rotating case, and 
they are therefore called "supramassive stars". These con- 
figurations are stabilized by rotation and would become 
unstable if slowed down below a critical rotation rate. 

So far we have considered stars chemical equilibrium at 
the center, i.e. fi n — Incidentally, for the EOS-class 
considered here, the resulting static configurations share 
a common outer surface in this case. However, global 
chemical equilibrium is generally not possible for config- 
urations with the two fluids rotating at different rates, 
which was shown by Andersson and Comer [9] and Prix 
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FIG. 6: Kepler configurations for EOS I. In figure (a) the relative rotation rate is TZ = 0.1, while in (b) it is 72. = 0.01. 
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FIG. 7: Gravitational mass M as a function of central baryon 
number density n c for EOS I with fixed central proton frac- 
tion of x p = 0.05. The four curves correspond to the non- 
rotating case, the co-rotating Kepler-configuration (TZ = 0), 
and two Kepler-configurations with relative rotation rates 
of TZ = 0.5 and TZ = —0.5 respectively. The circle indi- 
cates a static configuration with central chemical potentials of 
Mn = A*p = 0.3 rribC 2 . The box indicates the maximum-mass 
configuration in the static case. The dotted line represents 
the sequence of constant baryon mass connecting to the static 
maximum-mass configuration. Configurations above this line 
have no stable non-rotating counterpart and are called "supra- 
massive" stars. 



et al. [10]. In order to model more "realistic" configu- 
rations, in which the proton-fluid mimics a neutron-star 
"crust" (albeit without any solidity) by extending further 
outside than the neutrons, we can easily achieve this by 
choosing different central chemical potentials. For ex- 



ample, using EOS II and setting fj, n = 0.228 m^c 2 and 
/ip = 0.220 rrifcc 2 , we obtain the configuration shown in 
Fig. 8. For this figure we have chosen the rotation rate of 
the fastest known millisecond pulsar, which has a period 
of P r~j 1.56 ms. 
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FIG. 8: Configuration with protons rotating at the speed 
of the fastest known millisecond pulsar, Q p /2n = 641 Hz, 
and fl n /2-K — 645 Hz. The protons are extending further 
outside than the neutrons. The physical parameters are 
jU n = 0.228 rribC 2 and fi p = 0.220 mtc 2 , resulting in cen- 
tral baryon number density n c — 0.561 fm -3 , proton fraction 
x p = 0.09 and a gravitational mass of M = 1.39 Mq. 

A similar configuration with /i n = 0.28 m^c 2 and 
/i p = 0.3 mfcC 2 rotating at the Kepler-limit for a rela- 
tive rotation rate of TZ = 0.01 is displayed in Fig. 9. As 
can be seen by the cusp-formation, the Kepler-limit is 
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determined by the outer fluid, i.e. the protons in this 
case, despite the fact that they are rotating more slowly 
than the neutrons. This contrasts with the case depicted 
in Fig. 5, in which the faster fluid also happens to be the 
outer fluid, which is a particularity of this EOS-class and 
the choice of central chemical equilibrium fi n = fi p (cf. 
[10]). 




-10 10 



x [km] 

FIG. 9: Kepler configuration with TZ = 0.01, fi n = 0.28 m b c 2 
and fi p — 0.3 mjc 2 , corresponding to a central baryon number 
density of n c — 0.716 fm~ 3 , proton fraction x p — 0.09 and a 
gravitational mass of M = 1.57 Mq. The protons extend to 
the outer surface. The maximal rotation rates are found as 
Qn/2% = 924.5 Hz and n p /27r = 915.3 Hz. 



VII. CONCLUSIONS 

We have developed a theoretical framework and a nu- 
merical code for computing stationary, fully relativistic 
superfluid neutron star models. 

Using this code we have reconfirmed the existence 
of oblate-prolate shaped two-fluid configurations, previ- 
ously shown in [10, 14]. We have studied the dependency 
of the Kepler rate of a two-fluid star on the relative rota- 
tion rate TZ. We have compared this to the Kepler-rate 
predicted by a Newtonian slow-rotation approximation, 
which is found to typically overestimate the Kepler-rate 
by about 15%, but which can also sometimes underesti- 
mate it, as seen in the case of EOS I and TZ < —0.25 (cf. 
Fig. 5). 

The relative rotation rate can also have a large influ- 
ence (at fixed central density) on the mass-density rela- 
tion, as shown in Fig. 7. 

Another interesting aspect of this model is that we are 
not restricted to configurations in chemical equilibrium 
at the center. Choosing the central chemical potentials to 
be different allows one to emulate a neutron star "crust" 
(albeit a fluid one), as one fluid will now extend further 



outwards than the other, as seen in Fig. 8 and Fig. 9. 
One interesting observation from such configurations is 
that the Kepler-limit will be determined by the outer 
fluid (forming a cusp) , while this can actually be rotating 
slower then the inner fluid. 

We currently use a (quite general) EOS class of two- 
fluid polytropes, but this can be extended straightfor- 
wardly to more "realistic" nuclear-physics equations of 
state. In particular it might be interesting in the next 
step to use the first relativistic two-fluid EOS incorporat- 
ing entrainment by Comer and Joynt [13]. Furthermore 
it would be important to add the presence of a solid crust 
and to allow for differential rotation in the superfluid neu- 
trons (differential rotation in single-fluid stars has been 
implemented and used in LORENE already, cf. [29, 30]). 

The astrophysically most interesting future extension 
of this work would probably consist in studying the oscil- 
lation modes of such models, which would be directly re- 
lated to the emission of gravitational waves. In these non- 
stationary situations, however, dissipative mechanisms 
like viscosity and mutual friction would also start to play 
a role and should be included in the model. 



APPENDIX A: THE NEWTONIAN ANALYTIC 
SLOW-ROTATION SOLUTION 

A method for solving the stationary 2-fluid configura- 
tion in the Newtonian slow-rotation approximation was 
initially developed in [8], and was completed to include 
all EOS-interactions in [10] (in the following referred to 
as Paper I). Using this method, an analytic solution 
was found in Paper I for equations of state of the form 
£ = \n n n\ + TjKpTip + K np n n n p + /3 p n p A 2 . While this so- 
lution was very useful for studying the qualitative prop- 
erties of an interacting 2-fluid system, it is unfortunately 
not very suitable for comparison to the numerical solu- 
tion presented in this paper. The reason for this lies in 
the somewhat unphysical behavior of entrainment in this 
model. Namely, the entrainment numbers (16) are found 
as e p = 2/3 p /m, and e n = 2x p /3 p /(m(l — x p )), where 
x p = n p /n is the proton fraction, which is constant for 
this EOS (cf. Paper I). Therefore the entrainment num- 
bers are constant, independently of the densities, and so 
the entrainment effect would still be present in a 1-fluid 
region. This unphysical behavior does not pose a prob- 
lem in the slow-rotation approximation, which consists of 
an expansion around a static chemical-equilibrium con- 
figuration: the two fluids share a common surface in the 
unperturbed state, and the rotation will only induce in- 
finitesimal displacements of the fluids. In this framework 
there are therefore no finite 1-fluid regions. However, in 
a numerical code allowing for arbitrary rotations and de- 
viations from chemical equilibrium, such an entrainment 
model would be problematic. The EOS-class (68) used 
in this work is therefore preferable on both physical and 
numerical grounds. 

Fortunately, an analytic solution can also be found for 
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this physically preferable EOS using the slow-rotation ap- 
proach developed in Paper I. This solution is very valu- 
able for quantitative comparisons with our numerical re- 
sults presented in section VC. Here we derive this new 
analytic solution, skipping some of the more technical 
steps, which have been explained in more detail already 
in Paper I. 

Because of axisymmetry, the rotating solution only de- 
pends on the spherical coordinates r and 9, while the 
static configuration is assumed to be spherically symmet- 
ric. The 2-fluid slow-rotation approximation proceeds by 
expanding any local stellar quantity Q as follows: 

Q(r, e ; n x ) = o (0) (r) + n x cf Y (r, 9) n Y +o(n 4 ), (Ai) 

where here and in the following we automatically sum 
over repeated constituent indices X = n, p. We can sep- 
arate the variables r and 9 by expanding in Legendre 
Polynomials, i.e. Q^faO) = Qf y (r)P ; (cos (9), and 
it can be shown that only the components I = 0, 2 will be 
nonzero in the solution. The solution is therefore fully 
determined by two ordinary differential equations for the 
components $^ y (r) and <f>f y (r) of the perturbation of 
the gravitational potential. The information about the 
EOS enters via the following two "structure functions", 
defined as 



Sxy = 



d 2 s 



dnxdny 



d 2 e 



dn x dA 2 



(A2) 



where |o denotes the derivatives to be evaluated at the 
static configuration. For the EOS (68), we find 



5' / "'P ^np \ 

XY = F 

/V- I K np K n J 

where K, = n p K n — K^ p , and 

/3*(r) - K A M XY n { °\r), 
with the constant matrix M XY defined as 



(A3) 



(A4) 




(A5) 

We further introduce the "derived" structure functions, 

Ua = <Sabw b , and k = m A kA , (A6) 

which are constant for this EOS. The matrices E XY , de- 
fined as 

E XY (r) = Is*, (S b - xy 2/3» A XY ) , (A7) 

are now functions of r, contrary to the EOS treated in 
Paper I, in which they were constant. The constant aux- 
iliary matrices # 4 ' xy and A XY are defined as 



§ n,XY = / TO" \ 6Pi xY = [ 

\00/' I ™ p 



A XY = 



1 -1 
•1 1 



(A8) 



(A9) 



The static background solution only depends on S X y , 
and is identical to the one found in Paper I. Namely, 
in "natural units" defined by p^(0) = 1 and R=l, this 
static solution can be written as 



p<°>(r) 



sin 



(rVk) 



(A10) 



In these units it must be true that p(°)(l) =0, which 
leads to the condition k = m A kA = 7r 2 . This relation can 
be used to rescale the EOS parameters K n , n p and K np to 
natural units. The respective particle number densities 
nf\r) are expressible as 



nTir) 



^ A (0)r \ 



Substituting this into (A7), we can write 

E X A Y {r) = E XY -E XY P ^{r) 
in terms of the two constant matrices 



(All) 



(A12) 



E XY 



EV 



-Sab5 bxy 



3p(°) 



Sab^A 



XY 



(A13) 



and for E XY = m A E XY , we write in an analogous manner 
E XY = E XY -E XY p^{r), (A14) 

with 



& Y = h B s B * Y , 

E XY = ^Kk p A XY . 



(A15) 



We can now write the differential equations determining 
the solution for the given EOS, namely 

V ^ Y + tt 2 ^ y = C XY + r 2 ^ - rV )^ (A16) 
V 2 <f XY +n 2 <^ Y = -r 2 E XY +r 2 p^E XY , (A17) 

where the differential operator T>i is defined as 

d^_ + 2d__ 1(1 + 1) 

dr 2 r dr r 2 



V, 



(A18) 



We note that the only difference of this EOS to the one 
studied in Paper I concerns the entrainment (^(r). We 
can therefore formally recover the results from Paper I 
in the limit k a — > 0, which corresponds to E XY — ► and 
E XY — > 0. The constant matrix C XY is determined by 
the choice of stellar sequence, e.g. either characterized 
by fixed central densities (FCD) or fixed masses (FM). 
The solution to the above equations also determines the 
density distribution of the two-fluid star, namely via the 
relations 



"£o(»0 = SabC 



■BXY 



r 2 E XY 



kA^ 1 



n XY 
n A,2 



(r) = -^EY - k A ^ Y , 



(A19) 
(A20) 
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where the constants C A,XY are also determined by the 
choice of stellar sequence, and they satisfy the relation 
C XY = k A C AXY . The complete slow-rotation solution for 
the density distribution of the two fluids can be written 

as 

n A (r, 6) = n A °\r) + Q x (n% + n^P 2 (cos0)) tl Y . 

(A21) 

The general (regular) solution of equations (A16) and 
(A17) can be found explicitly as 



sequence: 



XY 6k A E XY ( sinr7r i r 2 n 2 

n A,0.FCD — i — 



rir 



1 ) + E XY r 2 



k A E XY ( 2 2 sinrvr 2 2 2 

— —. — (1 — r ir ) 1 — —r 7r ) cosr7r 

47T 4 V rir 3 ' 



E XY 



—rn smm , 



(A29) 



XY 



6 



E 



XY 



12tT 4 



$f» = A 



xy J 5/2(r-n-) E 



XY 



E 



XY 



127r 7 r 3 



C 



XY 



{ 3m sin nr + (3 - 2?r 2 r 2 ) cos rir } , ( A22) 



{ (45 + 27r 4 r 4 ) tit cos r7r 

+15(rV - 3) sinrvr} , (A23) 



k A E XY ( 2 5 J 5 / 2 (m) 



IT 

. o/, :i /-: vv r/- j 4 , , . 

+ - I — r — 3 I T7T COS T7T — [T TT — d) SlUTTT 



E XY r 2 



6 7r 5 r 3 I \ 5 



E XY 
H — r7r sin rir , 

7T Z 



(A30) 



in terms of the constant "structure matrices" E and E 
defined in Eqs. (A13) and (A15). 



where A XY and A XY are constants of integration, and 
J n {x) are the standard Bessel functions. One can ver- 
ify the asymptotic behavior $f Y <~ r 2 as r — > 0, which 
is required for regularity. In addition to the regularity 
requirements at the center, the solution must satisfy the 
following boundary condition at the surface (r = 1): 



$r'(i)+/(/+i)$r(i) = o. 



(A24) 



These boundary conditions result in the following rela- 
tions for the integration constants A XY and A XY : 

A^^2A XY = 12(n 2 - 2)^ + (1 - tt 2 )^ 

(A25) 

{3 + 2it 2 )E XY (A26) 



V2A 



XY 



+4n 2 C XY , 

—E XY — - 
tt 2 12tt 4 



Fixed central density (FCD) sequence 

The FCD-sequence is the most directly comparable to 
the numerical results discussed in this paper. This se- 
quence is defined by the condition n^Q(O) = 0, and in 
this case the remaining constant of integration can be 
determined as 



C 



XY 
FCD 



-3 1 



-E XY 
4 



and we also have the relation 
S A bC ¥CD 



■BXY _ k A p XY 
~~ 2~ L FCD 



(A27) 



(A28) 



Fixed-mass (FM) stellar sequence 

For completeness we also give the solution correspond- 
ing to a fixed-mass sequence, which might even be phys- 
ically more interesting. The difference to the FCD- 
solution only concerns the I = component, while rrj 2 
is the same in both cases. As discussed in Paper I, the 
FM-sequence is characterized by the conditions 

/ r 2 n XY tFM (r)dr = 0, (A31) 
Jo 

which lead to the following condition for the potential 

*£fm(1) = 0. (A32) 
This results in the integration constant 



C ™ = ( ^2 ~ 1 1 E 



while we can similarly determine SabC?^ from (A31). 
Inserting this into (A22) we get 



rpXY 

«mW = ^{r 2 -l + V2 



E 



XY 



+— —r ( 2tt 2 - 3 + (2?r 2 r 2 - 3) costtt 

127T 4 \ 

The I = density coefficient is therefore found by using 
(A19): 



Putting all the pieces together, we arrive at the following 
explicit solution for the density perturbations of the FCD 



n XY 
"A.O.FM 



kAE- 



XY 



, 30 + 3ir z — 5r 2 7r 2 smr7r 

57r 4 \ r i 



IOtt 
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+E XY I - 2 



r — 



3(1 



r ) — r7rsinr7r 



^36(1 - -^) - (3 - 2r 2 7r 2 )cosr^ 



7T" 



12tt 4 

— (1 + 3r 2 )— sinr7r ) , 
r / 



(A35) 



which completes the analytic solution in the FM case. 



Calculating the Kepler- 



We briefly review the method of calculating the Kepler- 
limit using the slow-rotation solution presented above. 
The result of this calculation was used for the Newtonian 
slow-rotation Kepler-limit presented in Fig. 5. As derived 
in Paper I, the Kepler-rate to order il 2 for each of the two 
fluids can be expressed as the solution of the equation 



where the zeroth-order expression is 

^o)=* (0) '(l), 
and the second-order correction terms reads as 



3 

k A 



71 XY 

n A,Q 



; n A,2 



For the EOS-class considered here, we find 



% } = -Gp(0). 



(A36) 



(A37) 



(A38) 



(A39) 



s<ti Y 



FCD 



9E* Y 



E 



•XY 



12n'- 



(6 - 7^ 2 ) 



2k A 

+ (5^ 2 - 24) , 



'XY 



(A40) 



6 fl 



FM 



E 



XY 



E 



iXY 



27 



4tt 6 



(216 - 39^ 2 + 2tt 4 ) - -t^EY . (A41) 



For each fluid A, we find the Kepler-limit J1k,a(^b) 
(where B ^ A) by solving the quadratic equation (A36). 
The Kepler-limit is then interpreted as the correspond- 
ing solution for the faster fluid, which in this case corre- 
sponds to the outer fluid, i.e. we have 



fiK,n(ft), for ft>0, 

n KtP (TZ), for TZ<0. 



(A42) 
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